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We present a tight-binding investigation of strained bilayer graphene within hnear elas- 
ticity theory, focusing on the different environments experienced by the A and B carbon 
atoms of the different sublattices. We find that the inequivalence of the A and B atoms is 
enhanced by the application of perpendicular strain e^z , which provides a physical mecha- 
nism for opening a band gap, most effectively obtained when pulling the two graphene layers 
apart. In addition, perpendicular strain introduces electron-hole asymmetry and can result 
in linear electronic dispersion near the K-point. When applying lateral strain to one layer 
and keeping the other layer fixed, we find the opening of an indirect band gap for small 
deformations. Our findings suggest experimental means for strain-engineered band gaps in 
bilayer graphene. 



I. INTRODUCTION 



The synthesis of graphene — a single layer of carbon atoms arranged into a honeycomb structure 
— and measurements of its electronic properties in 2004 by Novoselov et al. [T] has sparked off the 
development of a whole new research field, continuing to expand rapidly today in both experimental 
and theoretical directions. The high rate at which the graphene literature is growing is reflected by 
the regular appearance of reviews on graphene in recent years; for a selection of relevant accounts 
we refer to Refs. [2H8]. 

The original excitement about graphene not only came from its two-dimensionality [H |9] — flat 
two-dimensional (2D) crystals had not been successfully fabricated before and were even predicted 
to be unstable — but also from its electronic properties [H [9l-[TT] . Apart from novel fundamental 
physics, graphene boasts superior material properties, including high thermal conductivity, current 
density, carrier mobility, carrier mean free path, strongness, stiffness, elasticity and impermeability. 

Not only single graphene sheets but also stacks of a few graphene layers [1] where the layers are 
coupled by van der Waals interactions as in graphite can be isolated. This has lead to investigations 
of multilayer graphene, and in particular of bilayer graphene (Fig. [T| top). Interestingly, bilayer 



2 



ff = 2 




E 
t 




K -^k 



K -^k 



FIG. 1: Graphene bilayer, with A carbon atoms in black and B carbon atoms in gray (top). Electronic 
spectrum near the K-point in absence (bottom left) and presence (bottom right) of a bias voltage W between 
the layers. 

graphene displays (almost-)parabolic electronic dispersion at the K-points (Fig. [T| bottom left), 
making electrons behave differently |12H15j as compared to the single-layer case. Bilayer graphene 
offers the possibility of applying a bias voltage W between the two layers, allowing to tune the 
band structure. In particular, the inequivalency of the two graphene layers then gives rise to a 
"Mexican-hat-like" band structure featuring a band gap I15H19I of magnitude E„ = l*^^!*-!- 
(Fig. [T| bottom right), with t_L ~ 0.377 eV the interplane hopping parameter (see below). A 
tunable energy gap is important for possible electronic devices. 

Another means of influencing (mono- or bilayer) graphene's electronic structure, currently re- 
ceiving a lot of theoretical attention [20H28j . is provided by mechanical deformations. It has e.g. 
been shown that it is possible to conceive inhomogeneous strains in single-layer graphene such that 
they act as a high uniform magnetic field, therefore resulting in strain-induced Landau levels and a 
zero- field quantum Hall effect |20J . In Ref. |21| , Pereira et al. elaborated a tight-binding description 
for uniaxially strained single-layer graphene and predicted that a band gap can form upon defor- 
mations — along preferred directions — beyond 20%. Recent investigations on strained bilayer 
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graphene include a standard tight-binding treatment of uniaxial strain |22j . a description of elas- 
tic deformations and electron-phonon coupling in bilayer graphene by means of pseudo-magnetic 
gauge fields the effect of strain on the Landau level spectrum and the quantum Hall effect 
\24:\ [25] and the effect of strain in combination with an external electric field [26[ [27] . 

In the present paper, we employ a nearest- neighbor tight-binding description and linear elasticity 
theory to show that a perpendicular strain component modifies the on-site energies of the two 
carbon sublattices in bilayer graphene, which can open a band gap at the K-point. The band 
gap is of a different nature than in the case of the "Mexican-hat-like" electronic dispersions. In 
addition, we consider the case of two graphene layers subjected to different (but uniform) strains, 
a scenario proposed recently and predicted — by means of ab initio calculations — to also lead 
to the opening of a band gap [28] . Uniform deformations parallel to the sheets are studied as well 
and are compared to the monolayer case |21) . 



II. STRAINED BILAYER-GRAPHENE 



We first consider the unstrained AB (Bernal) stacking variant of bilayer graphene: two graphene 
layers (labeled 1 and 2) at c = 3.44 A apart with the A atoms of layer 2 sitting directly on top of 
the A atoms of layer 1 (Fig. [T| top). The B atoms of layer 1 and the B atoms of layer 2 have no 
direct neighbor in the opposite layer. 

The band structure of bilayer graphene can be described within the tight-binding formalism. 
The nearest-neighbor tight-binding hamiltonian assumes one free 2pz electron provided by each 
carbon atom and reads 

^ Xc Xc ' \^ Xi+di Xi+di X2-di X2-di I 

" \Xi X2 I 

/ 3 3 _\ 



•Xi X2 



1=1 



+ t±y^ fct , + ct ^C^ 1 . (1) 
Z-^ I Xi Xi+c Xi+c Xij ^ ^ 

Xi 

The index a stands for the layer (1 or 2), the vectors are the lattice sites of the A atoms 
of layer a. Every A atom in layer 1 is surrounded by three B atoms at relative position vectors 
di = ^acx + ^aiy, d2 = \aex — ^cL^y and djj = — ae^, where the bond length a has a value of 
1.42 A, and the A atoms in layer 2 have neighboring B atoms at relative positions —di, / = 1, 2, 3. 
The operators and create and annihilate an electron at site X, respectively. The vector 
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(c = 3.35 A) connects two nearest-neighbor atoms in different graphene sheets. Values 
for the intra- and inter-plane hopping energies t and are obtained from fitting the tight-binding 
model to experimental data for graphite. Here we use the values quoted in Ref. [29J: i = 3.12 eV 
and = 0.377 eV. The on-site energies Va and differ slightly due to the different environments 
of A and B atoms. The difference A = |Va — Vb| ~ 0.009 eV [29j is about two orders of magnitude 
smaller than the hopping parameter and is usually considered only in models going beyond the 
nearest-neighbor tight-binding hamiltonian ([T]) where also A1-B2 and B1-B2 hoppings are taken 
into account (see e.g. Ref. (HJ). We therefore put Va^Vb^V ■ The value of the on-site energies 
V will turn out to be of critical importance when considering perpendicular strain; we will return 
to it later. For a comprehensive review on a complete tight-binding description of (unstrained) 
bilayer graphene, and in particular for a discussion of features in the electronic structure resulting 
from asymmetry of the diagonal (differences in on-site energies, Va 7^ Vb), we refer to Ref. j30j . 

The direct-space hamiltonian ([T]) can be converted into a reciprocal-space hamiltonian by in- 
troducing four-component spinors 



= ( 4(gl, 4(gl, 4(0, 4(gl ) , 



(2a) 



(2b) 



(2c) 



Here, the operators Ci{q) {i = Ai,Bi, A2,B2) are the discrete (lattice) Fourier transforms of c 



X,: 



X, 



■iq-Xi 



(3) 



where the ^-vectors of the first Brillouin zone are defined so that the properties Ylq^^'^ ~ 
and X^j? e^"^'^' = NS^^^ hold. The hamiltonian matrix Sj{q) reads 



( V c{q) i± \ 

C{q)* V 

t± V Ciq}* 

V C{q) V J 



(4) 
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with 



do) 



1=1 



(5) 



The band structure {E{q)} is obtained by solving the secular equation det(^S){q) — E^qjl^j = 0, 
with /4 the 4x4 unit matrix, in the first Brillouin zone. At the K and K points, two electron 
energy bands touch each other at the Fermi level, making the material a semi-metal (zero band 
gap). 

Since the electronic properties of (bilayer) graphene are determined by the band structure near 
the point — qk = |f e*c + g^^ e^ and q^i = |f — ^^^ e^ — , it is convenient to make a Taylor 
expansion of f){q) around q^o- With q = q]^ + k and retaining only lowest-order terms in k, the 
hamiltonian matrix near becomes 
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(6) 



Here, (j) is defined via kx + iky = ke^'^'. For q = q^ + k, the hamiltonian matrix {k) is the 
complex conjugate of S^^{k). Note that the quantity ^ can be rewritten as hv-p, with vp the Fermi 
velocity. Diagonalisation of ^ results in the dispersion shown in the left bottom of Fig. [T] 

As mentioned in the introduction, the application of a bias voltage W between layers 1 and 2 
(replacing the elements i^f/ ^ {k) and i^^g ' {k) by eW/2 and i^^j ^ (^) and S)^^ ' (k) by —eW/2) results 
in a finite band gap p^HT9] . However, this band gap lies not at g = ^j^o but at a 5- vector slightly 

away from q-^o (Fig. [T| bottom right), and its theoretical maximum value is \\mw ^^o Eg = t±. 

In the following, we will show that applying a strain to the bilayer graphene lattice can also lead 
to a band gap, which is not bounded by t±. 

In the presence of a displacement field u{X), the position of an atom formerly at X is X + u{X). 
The effect of a small deformation of the lattice can be written as a correction 6H to the original 
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hamiltonian H, with 



6H 



'^^A E E + ( E + E 

+ E'^*± 



-di + 



(7) 



Within the nearest-neighbor tight-binding approximation, the changes in on-site and hopping en- 
ergy parameters Va, Vb, t and t± are related to changes in nearest-neighbor interatomic distances 
(bond lengths). We stress that it is therefore important to distinguish between A and B sites since 
they have different environments in the bilayer. As pointed out before, an A site has three in-plane 
nearest-neighbor B sites and one neighboring A site in the opposite layer at a distance c + 5c; 
a B site has only the three surrounding in-plane A sites as nearest neighbors (see Fig. [T| top). 
Denoting the three (not necessarily equal) changed bond lengths between neighboring in-plane A 
and B atoms by a; = a + 5ai (/ = 1,2,3), we have for the corrections to the on-site energies to 
linear order in the deformations 
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6ai + 



dV 
dc 



5c, 



da 



5ai 



5^3 = E 

1=1 

for each of the two layers. For the hopping parameters we have 

dt ^ 
oti = T- 5ai, 

5t± = ^i:^ 5c. 



(8a) 
(8b) 



(9a) 
(9b) 



In the following we will drop the attributes \a and \c- The link between the corrections 5Va, 5Vb, 5ti 
and 5t± and the displacement field u{X) then comes from considering the bond length corrections 
5a and 5c. Details of the calculations and approximations involved are given in Appendix [Aj the 
resulting correction 5S){q) to the hamiltonian matrix Sj{q) reads 



3a aV 








3a eV 
2 da 
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(10) 
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with 

3 

1=1 

'^5^^ U2 ^°^^^^^y^ + ^ + - cos{—aqy)eyy + V3^sln(— agj,)e^y 

(11) 

Here, the quantities £ij = x,y,z) are elements of the strain tensor [e], the general definition 
of which involve derivatives of the displacement field components to the coordinates. For the 
uniform displacements we mostly consider in this work the elements £ij can be related to relative 
increments/decrements of the bond lengths occuring in the lattice (see Appendix|A|: u{X) = [sjX. 

The leading correction 6S)^ to the hamiltonian matrix Sj^{k) at the K-point is /c-independent 
and reads 



•3a Si/ I _ OjV ^ 3aiW( 



^ dc 
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(12) 



For q = + k, the hamiltonian matrix correction 5S^^ is the complex conjugate of 5S^^. The 
K-point hamiltonian correction for a strained single graphene layer, of the form 



(13) 



9l{^xx ~l~ ^j/y) 92{^xx ^yy ~l" '^'^^xy) 

has been derived before in the context of electron- phonon coupling in carbon nanotubes |31] . 
where the term proportional to gi is a deformation potential, a concept going back to Bardeen 
and Shockley |32j and the term proportional to §2 corresponds to a bond-length change. However, 



to our knowledge, the derivation of the strained bilayer hamiltonian [Eqs. (10) - (11)] as given 
in Appendix [A] — with particular emphasis on the asymmetry on the diagonal — has not been 
reported before. 

III. PERPENDICULAR UNIFORM STRAIN 

We first consider the bilayer-specific possibility of perpendicular strain, Szz, associated with a 
change from the inter- layer distance c to c' = c{l + ezz)- Putting the in-plane strain tensor elements 



8 



/ = 0.5 




-0.3 -0.2 -0.1 0.1 0.2 0.3 -0.3 -0.2 -0.1 0.1 0.2 0.3 -0.3 -0.2 -0.1 0.1 0.2 0.3 -0.3 -0.2 -0.1 0.1 0.2 0.3 -0.3 -0.2 -0.1 0.1 0.2 0.3 



fc(A- 



k (A- 



k (A- 



k (A- 



k (A- 



FIG. 2: Low-fc band structure E(k) for perpendicularly uniformly strained bilayer graphene for various 



values of \c^ez 



to zero, the hamiltonian Sj^{k) + 6Sj^ takes on the form 



hvpke 



dc 

i<f) 



dc 












V 



g^ozz hvpke 
hypke^'t' V 



-i(f> 



(14) 



V 




The phase factors e^*6 present in Eq. ^ have been eliminated by means of a redefinition of the 
operators for B sites [Eqs. ( |A24a[ ) - ( |A24b[ )]. 



The symmetry-breaking along the diagonal exhibited by the hamiltionian matrix ( 14 ) is unusual 



since it distinguishes not between layers (as e.g. a bias voltage between the two layers would do) 



but between A and B sublattices. Interestingly, the hamiltonian (14) displays the possibility of the 
opening of a band gap at /c = 0: if 



dV 
oc 



dt± 

> t±+ c—-ezz, 
oc 



(15) 



the degeneracy of the bands at A; = is lifted, as illustrated in Fig.[2| (Criterion (15), an analytical 

is small compared to t_L.) Recently, Mucha- 



result, holds when t± + c^^ezz > 0, i.e. when c 
Kruczyhski et al. [30j examined the asymmetry of the (unstrained) graphene bilayer hamiltonian's 
diagonal. The possibility of band gap openings and electron-hole asymmetry, as encountered here 
in Fig. [2| due to different on-site energies, was realized. Here, we show that strain enhances 
the diagonal's asymmetry [Eq. ( |10[ )] and that elastic deformations therefore provide a physical 
mechanism for the band structure modifications discussed in Ref. |30j . 

To check whether the possibility of a strain-induced band gap is experimentally relevant, a more 



dc ^ 



quantitative investigation of criterion ( 15 ) is in order. First, the variation of the hopping parameter 



t± with interlayer distance c can be estimated using Harrison's relation 



dt± 

dc 



(16) 
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FIG. 3: (Color online) Plots of the variation of ^ 2^ (blue, full line, expansion) and ^ + 2^ 

(violet, dashed line, contraction) with Ezz- The two horizontal lines represent the values of |^| resulting 
from using the two extreme values for V considered in the present work — V — —2 eV (olive, dashed-dotted 
line) and V — —0.25 eV (green, dashed line). 



which fohows from an assumed inverse-square dependence of t± on c |33j . The inequahty (15) then 
becomes 



dV 



dc 



(17) 



For expansion (ezz > 0), this can be rewritten as 

dV ^ 



dc 

while for contraction (ezz < 0), one obtains 

dV 



1 



(18) 



dc 



t± 



+ 2 



(19) 



The intra- and inter-plane hopping parameters have values of t = 3.12 eV and t± = 0.377 eV, 
respectively. The interplane distance is c = 3.35 A, so that ^ = 0.113 eV/A. Using these values, 

we have plotted the quantities ^ (^-f" 2^ and ^ (|F~7 ~'~ ^) -^^S- j^j Interestingly, it follows 

that the band gap formation criterion is reached more easily (smaller \£zz\) in the case of expansion. 

To obtain a physically meaningful estimate for the quantity |^| we proceed as follows. First, 
we recall the tight-binding definition of V: 



V = £2p^ = J dfiP*ir)U{r)^{r}. 



(20) 



Here, ip^f) is the 2pz electron wave function for a carbon atom, and U{r) is the periodic potential 
of the lattice. Suprisingly, while the tight-binding formalism is a standard method for modelling 
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band structures — especially for carbon structures — , a numerical value for V = e2p^ has not yet 
been established in the literature. The plausible reason is that when the difference between Va 
and Vb is neglected (see Sect. |ll]), the presence of V on the diagonal of the hamiltonian matrix 
merely results in an overall shift of the energy bands, making its actual value redundant. In Ref. 
|34j . Reich et al. present a detailed comparison of tight-binding and ab initio energy dispersion 
calculations for graphene and provide fits of the tight-binding parameters to the band structures 
obtained by density-functional theory. Two variants were considered: one fitting the full MTKM- 
path in /c-space and one where only A:-vectors yielding optical transitions with energies less than 
4 eV were taken into account, resulting in e2p^ = —0.28 eV and —2.03 eV, respectively. We will 
therefore consider the range —0.25 eV > e2p^ > —2 eV. The simplest way of modelling U (r) is to 
put positively charged ions (charge le > 0) at r = and r = c'e^: 

?7(r) PC - ,^ \^ , = -- - , ^ (21) 

1^ \r-de^\ r ^^2 _ 2dr cos 6 + c'2 

where spherical coordinates have been introduced. For the carbon 2pz electron wave function, we 
follow the common practice of taking the hydrogen-like 2pz orbital: 

'4}{r) oc re ^"o cos0, (22) 

where ao = 0.529 A is the Bohr radius. We then get for V the expression 

V = c[ r'^dr [ sin OdOr'^ e~ ^ cos"^ 9 (-- - , ^ |, (23) 

Jo Jo \ r Vr2-2c'rcos^ + c'2 7 



where the proportionality factors left unspecified in Eqs. (21) and (22), together with the factor 
27r coming from the azimuthal integration, have been collected into the factor C. The integrals 

^ / N r in sin6'cos2 6l . ^, 

SnAr)= dO (24) 
Jo V ?^ — 2nc r cos + [ndy 



entering Eq. (23) can be solved analytically: 



So{r) = |:, (25a) 



2 5(rac')^+2r^ 



if < r < nd 



T% -^^ — T \i r > nc 
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For c' = c = 3.35 A, the integral 

D{c) = / drr^e «o [S^ir) + Sic{r)) (26) 
Jo 
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FIG. 4; Plots of the variation of V{c') with c' = c(l + e^^) for V{c) = -0.25 eV (blue, fuU hne) and 
V{c) = —2 eV (violet, dashed line). 

has the value D{c) = 0.539 A^. From V = CD{c) we can then obtain the "calibration value" 
C = V/D{c). For the values V = -0.25 eV and V = -2 eV, we obtain C = -0.463 eV/A"^ 
and C = —3.709 eV/A~^, respectively. The dependence of V on c' reads V{c') = CD{c'), it is 
visualised in Fig. [4} Clearly, a linear approximation in the range —0.25 < < 0.25 is justifiable, 
and the value of ^ can be easily determined numerically. We find ^ = 0.0287 eV/A and 0.230 
eV/A for V = —0.25 eV and V = —2 eV, respectively; these values are marked by horizontal lines 
in Fig. [3| 

Within the present model, it follows that in the case of a large tight-binding on-site energy 
parameter (\V\ = 2 eV), a band gap opens for positive strains larger than Szz ~ 0.25 (see Fig. 
[3]), corresponding to interlayer distances larger than c' ~ 4.19 A. For Szz ~ 30%, the band gap's 
magnitude is about 125 meV (Fig. [2j / = 1.5). 

We recall that formally, criterion ( jls] ) is valid when the right-hand side is positive, i.e. when 
Ezz < 0.5, hence the upper limit of £zz 

= 0.5 in Fig. [3j Going beyond strains of 50% (both positive 
and negative) would be well outside the validity of the assumption of small displacements u{X), on 



which the hamiltonian matrix (10) relies (see Appendix [A|) . For e ~ 0.25, neglected contributions 
are of the order 0.063 which is still acceptable. Similarly, we point out that the unphysical 
behavior of an ever increasing band gap with increasing Ezz must become invalid when linear 
elasticity, i.e. the assumption of small bond length changes [Eq. ( |A2p ], fails. 

Based on the foregoing elaborations, stating that the opening of a strain-induced band gap by 
pulling apart the two graphene layers may be experimentally observed is a fair conclusion. Our 
main purpose here is not to provide accurate predictions but rather to point out the consequences 
of the symmetry-breaking along the diagonal in the hamiltonian of strained bilayer graphene. 
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Accurate density- functional theory calculations could provide better numerical estimates, but most 
importantly, experiments should be undertaken to investigate the effect on the band structure upon 
pushing together or pulling apart the graphene layers. 

IV. SYMMETRIC UNIFORM xy-STRAIN 

We next consider pulling or pushing the bilayer in a direction parallel to the graphene sheets 
i^zz = 0). The tight-binding hamiltonian then reads 



f v + f^ie.x + syy) a<t) + sa<r) 

tx 




t_L \ 



(CD + KiQ) V+f^ie^. + Syy) J 



(27) 



The 4 solutions of the corresponding secular equation are 



t\ + 2\C{q) + bC{q)\^ ± t^Jtl + 4|C(g1 + Km 



(28) 



There can only be a band gap when (^{(pj + bQi^q) differs from zero for all g-vectors. The same 
condition arises from the monolayer strain problem, where the 2 solutions of the secular equation 
read E±{q) = ±\C{q) + KioJl (leading to the Dirac cone at the K*^')-point for 5Ciq) = 0). Pereira 
et al. PT] have investigated the behavior of \C{q) + SC{q)\ in detail. Using the isotropy of a 2D 
hexagonal lattice, the strain tensor can be written as 



cos^ 9 — a sin^ 6 (1 + cr) cos 6 sin 9 
(1 + 0") cos^sin^ sin^ 9 — a cos^ 9 



(29) 



with 9 the angle between the tension T and the y-axis (T = Tcos(7r/2+^)e2.+Tsin(7r/2+0)ey) [35], 
and a the Poisson ratio, which takes the graphite value of 0.165 which we choose in the remainder 
|36j . The conclusions made by Pereira et al. are that (i) the minimal strain required for opening a 
gap is about 23%, that (ii) tension along the zig-zag direction (9 = 0, 7r/3, 27r/3) is optimal for the 
opening of a band gap, and that (iii) tension along the armchair direction {9 = tt/2, 5tt/6, llvr/G) 



never results in a band gap. From Eq. (28) it follows that the inter-plane coupling t± does not 



play any role and that the same conclusions are valid for the bilayer. 



13 



V. ASYMMETRIC UNIFORM xy-STRAIN 

Finally, we consider the possibility of applying different strains (parallel to the bilayer) to the 
two layers a = 1 and a = 2. In Ref. [2H], ab initio calculations of a graphene bilayer's band 
structure were performed for the situation where one of the two layers is subjected to positive 
strain along the zig-zag or armchair direction (without the elastic response in the perpendicular 
direction). Here, we consider a more general situation. We choose to let layer 1 intact ([e] = [0]), 
while layer 2 is pushed or pulled ([e] ^ [0]). Any pushing/pulling direction is allowed; the elastic 



response is taken into account by using expression (29) for the strain tensor of layer 2. 

The deformation of layer 2 introduces a mismatch between the two honeycomb lattices. For 
infinitesimally small deformations e, the crystallographically correct unit cell of the bilayer — the 
"least common multiple" of the unit cells of the two layers — becomes infinitely large. For practi- 
cally feasible descriptions of the asymetrically distorted bilayer one is forced to choose deformations 
that lead to not-too-large unit cells, resulting in a discrete set of strain values e. For example, for 
their ab initio calculations, Choi et al. [28j considered minimal zig-zag strains of about 2%, re- 
quiring a supercell of 51 unit cells in layer 1 and 50 in layer 2. A tight-binding description faces 
the same mismatching problem. Apart from the necessity to set up a larger hamiltonian matrix, 
the inter-layer hopping parameters would have to be carefully reconsidered. Indeed, in the case of 
asymmetric strain, the A2 atom originally directly above a particular Ai atom now has a different 
position, while a Bi atom may now have a direct (or almost direct) neighbor above. In a way, the 
structure becomes a mix of AB and AA stacking |28j. A similar issue was addressed recently in 
Ref. [37J: spatial modulations in the bilayer interlayer hopping arising due to elastic shearing or 
twisting were shown to lead to a non-Abelian gauge potential in the description of the low-energy 
electronic spectrum. 

Designing the full tight-binding matrix with correctly interpolated inter-layer hopping parame- 
ters is a formidable task, the outcome of which would still only be a limited set of finite accessible 
strains e. As a compromise, we therefore consider the following small-A: hamiltonian: 



( V hv-pke^'l' t± \ 

hv^ke-^'t' V 

t± V +^^(e:,-^+eyy) hvpke-^'l' + §^{-e:,^ + Eyy - 2ie:cy) 

V hvT.ke^t' -i^§L{-e^,+eyy + 2ie^y) v + f ^{e^^ + Syy) 

(30) 
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where the parameter t± has to be interpreted as representing the average hopping between the two 
layers. The apphcation of asymmetric plain affects both the diagonal and the off-diagonal elements; 
the resulting band structure is therefore non-trivial and worth investigating. Diagonalising the 



hamiltonian (30) results in the following secular equation: 



{E{kf - IC(^)P) Uf - E{k)Y - m) + 5CP1 + tlE{k){F - E{k)) = 0, 



(31) 



with 



F 
5C 



3a dV 



2 da 



3a dV 
,3a dt 

~^^~da ^^^^ + ^yy + ^'''^xy 



e{l-a) 



3a dt 
Ida 



e[i{l - a) cos(26') + (1 + a) sin(26l)] , 



(32a) 
(32b) 
(32c) 



tan' 



-I ky 



does, in general, not cancel out, and that the full 2D 



Note that the argument ^ 
dependence of E{kx,ky) has to be considered. Note also that upon replacing e by —e in Eqs. (31) 



- (32c), the secular equation remains invariant if E changes sign and the phase (j) is shifted by 
vr. Hence, the substitution e — > —e only leads to a band inversion and it suffices, as far as the 
opening of a band gap concerns, to consider e > 0. 

The change in band structure comes from the interplay between the values of a^, a^ and t±. 



To obtain an estimate for ^ we use the equivalent of relation (16): 



m 

da 



-4.394 eV/A. 



(33) 



Note that Harrison's relation [Eqs. (16) and (33)] should be taken as a rule of thumb rather than 
as an exact result [33j. Other (experimental or theoretical) values than 2 for rj = — circulate 



in the literature (e.g. rj = 3.6 |38], r] = 3 f39] or 77 = 1.1 [iQ]). As it is our aim to make qualitative 
conclusions, we have used rj = 2. 



For we proceed as for ^ in Sect. HI and put a charge at r = and f = a'e, with e any 
unit vector perpendicular to the z-axis (a convenient choice is e = e^) to mimic the periodic lattice 
potential : 



U{r) oc 



1 

|r| 



1 



1 



a'e\ 



\Jt^ — 2a' r sin 9 cos </> -|- a'^ 



Using the definition (20) for V{a'), we numerically calculate the derivative 



0.0439 eV/A and 0.351 eV/A for the cases V{a' 



da \a'=a 

0.25 eV and V{a' = a 



(34) 

and obtain 
= -2 eV, 



respectively. Inserting the value 



da 



0.351 eV/A, we observe the opening of a band gap for 
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arbitrarily small values of e. In Fig. [5j low-A; band structures are shown for a few choices of {e, 0). 
The {kx,ky) electronic dispersions always consist of two facing doubly-peaked surfaces. This can 
result in an indirect band gap [Fig. [5]^a)], but also in indirect band crossings [Figs.[5];b) and (c)]. 
An inspection of the ranges < e < 0.25 and < < n allows to conclude that (i) above a 
critical value for the strain, only indirect band crossings are observed and (ii) the angle 9 has 
little or no influence on the presence/absence of a band gap. The former conclusion agrees with 
the observation of Choi et al. [28] of a decrease of the band gap for zig-zag strains larger than 
e ~ 9%. As for the dependence of the band gap on 0, Choi et al. [28j found that in the case of 
armchair strains, no band gap appears at all. We recall that we allow for an elastic restoring force 
perpendicular to the pulling/pushing direction which results in a dependence between Exx^ ^xy and 
Eyy [Eq. (29)]. In our model, this isotropy makes the band gap development independent of 6. 

A true (but indirect) band gap is only observed for small strains e, as e.g. in Fig. [§a) (e = 0.05, 
6 = 0, AEg ~ 10 meV). The maximal band gap turns out to depend on the magnitude of '^^ 



da ' 

for transversal strain where the value of ^ is critical. For 

oc 

the value ^ = 0.0439 eV/A, corresponding to V{a) = —0.25 eV, the band gap for e = 0.05 and 



not unlike the situation in Sect. 

av 

da 



III 



: reduces to AEg ~ 1 meV). For the artificially large value of ^ = 1 eV, we find (Fig. 
6|) Ai?g = 25 meV. At this point we remark that one can identify with the deformation 



potential gi entering the strained monolayer hamiltonian (13), whose value for graphite is about 
16 eV |41j . For ^ one therefore has approximately 7.5 eV A, which suggests that larger values of 
^ may indeed be more appropriate. Obviously, experiments and/or ab initio calculations should 
be carried out to clarify the parameters' values. 



In conclusion, hamiltionian (30), a possible model for describing the effect of asymmetric lateral 
strains in bilayer graphene, results in indirect band gaps for small strains {e ~ 5%). For larger 
strains (e > 10%), overlapping indirect bands result in a gapless spectrum. No critical directional 
strain dependence is observed, which comes from the anisotropy of the 2D hexagonal lattice. 



VI. CONCLUSIONS 



Applying strain to a graphene bilayer alters its electronic structure. Using the simplest tight- 
binding model, featuring only nearest-neighbor intra- and inter-plane hopping (t and t±), we have 
shown that for the AB (Bernal) stacking of two graphene sheets, several band gap scenarios are 
possible. 

A first possibility is to push (pull) the bilayer's two graphene sheets towards (away from) each 
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(c) 



FIG. 5: Band structure of asymmetrically in-plane deformed bilayer graphene near the K-point, for ^ = 
0.351 eV/A: (a) e = 0.05, 6* = 0, (b) e = 0.1, 6 = tt/S, (c) e = 0.15, 9 = tt/I. 

other. As a consequence of the different neighborhoods experienced by A and B carbon atoms, the 
change in on-site tight-binding energies results in a symmetry breaking which, for large enough 
strains, opens a band gap. We find that pulling and pushing are inequivalent; the former is 
more effective in producing a gap. Predictions of the minimally required strain to open a band 



17 




FIG. 6: Band structure of asymetrically in-plane deformed bilayer graphene near the K-point, for ^ = 1 
eV/A, e = 0.05 and 6* = 0. 

gap critically depend on values for the changes in on-site and inter-layer hopping energies with 
inter-plane distance and ^^)- These quantities are not readily available from the literature. 
We therefore rely on rules of thumb to make acceptable estimates, and find that for pulling, a 
band gap opens from Ezz = 0.25 onwards. For Ezz ~ 30%, the band gap is about 125 meV. 
For pushing, the critical strain is well beyond 50%. Note that a strain of 50% would require a 
pressure of p = 0.5c44 2 GPa (taking for the elastic constant C44 the value of graphite, C44 = 4.18 
GPa [42j), which is experimentally feasible. While bringing the two layers close together can 
be realized by applying high pressure, pulling the graphene sheets away from each other is an 
experimental challenge. A possible indirect way to do so would be to intercalate the bilayer; 
recently, for example, Li atoms have successfully been intercalated between two graphene sheets 
|43| . Interestingly, we find that at the critical strain, the electronic dispersion near the K-point 
consists of a triple degeneracy of two crossing linear bands (Dirac cone) and one parabolic band. 
This result shows that under certain symmetry-breaking strain conditions, bilayer graphene can 
exhibit Dirac fermions. Normally, Dirac fermions only occur in graphene multilayers consisting of 
an odd number of graphene sheets [H] . The asymmetry of the hamiltonian matrix diagonal lies at 
the basis of the band structure modifications observed in Fig. [2] It was already recognized before 
|30j that intra-layer on-site energy differences can lead to band gaps of a different type than the 
"Mexican-hat-like" band gap band structures |15H19j associated with inter-layer energy differences 



18 



e.g. induced by the application of a bias. In the present work, we have shown that uniform 
perpendicular strain provides a physical means for enhancing the intra-layer energy differences. 

A second possibility is to push/pull parallel to the bilayer (symmetric uniform strain). It turns 
out that the criterion for opening a band gap is the same as for monolayer graphene; the inter-plane 
coupling t± does not play any role. The conclusions made by Pereira et al. \21\ for the monolayer 
can be transferred to the bilayer; a strain along the zig-zag direction larger than ~ 23% results in 
a band gap, while pushing/pulling along the armchair direction never leads to a gap. 

Finally, we considered asymmetric uniform strain: the two graphene sheets experience different 
strains parallel to the bilayer. Assuming a tight-binding hamiltonian matrix where t± represents 
an average inter-plane hopping, we obtain the interesting result that a small finite strain, in any 
direction, immediately opens a band gap. The maximal band gap depends, similarly to the case 
of transversal strain, on the interplay of the quantities ^ and ^ , reliable estimates for which are 
hard to derive. The observed band gaps are relatively small — ~ 10 meV. The advantage, 
however, is that there is no strain barrier. Rather, beyond a certain strain (typically 10%), bands 
start to indirectly overlap and destroy the gap. As opposed to the case of uniform symmetric 
strain, no dependence on the strain direction was obtained — a consequence of the isotropy of a 
hexagonal lattice and of taking the elastic response in the perpendicular direction into account. 
Pushing and pulling are equivalent. We point out that the observed band gaps are indirect, which 
is relevant for possible opto-electronic applications. Our results are in qualitative agreement with 
recent ab initio calculations of the electronic structure of similarly asymmetrically strained bilayer 
graphene [28] . 

Although having provided estimates for strains e required for the opening of a band gap and 
associated band gap magnitudes AE'g, we wish to emphasize the qualitative aspect of our results. 
Irrespective of the precise values of the tight-binding parameters and derived quantities, it is the 
particular structure of the graphene bilayer that, when deformed, allows for gapped electronic 
structures. In our opinion, the various ways shown here in which a band gap can be induced in 
bilayer graphene by deformations should stimulate experimental investigations on the possibility 
of strain-engineering bilayer graphene's electronic properties. Recalling graphene's high strength, 
this should be feasible. In addition, we suggest that the theoretical models presented here be 
reconsidered by means of precise ab initio calculations. 
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Appendix A: Strain 



As described in the main text, deformations of the bilayer graphene carbon network enter the 
tight-binding formalism via bond length changes [Eqs. (8a) - (9b)]. In this Appendix, we elaborate 



the expressions for 6ai and 6c and the strained graphene bilayer tight-binding hamiltonian. 



1. Continuous description. Long-wavelength limit 



In general, one has 



a + 5ai 



X + di + u{X + di) -{X + u{X)) 



+ 2di ■ {u{X + di) - u{X)) + \u{X + di) - u{X)\ . 



(Al) 



Assuming small bond length changes 6ai, we retain only first-order terms: 



so that 



a + 6ai^a(^l + ^di ■ {u{X + di) - u{X))^ , 



6ai^^di-{u{X + di)-u(X)). 



(A2) 



(A3) 



In the case of small bond length changes, the differences between displacements at neighboring 
sites \u{X + di) — u{X) | must be small as well. In other words, the displacement field n is a slowly 
varying (vector) function of X, implying that a linearisation of n at X approximates u sX X + di 
well enough. In this so-called long-wavelength limit we then have 



di ■ [u{X + di) - u{X)) « {di ■ V){di ■ u) 



(A4) 



The operator V = acts on u{X) with X taken as the continuous spatial variable 

X = {x,y,z). Writing di^, a = 1,2 for the x- and y-components of di (di has no z-component, see 
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top of Fig. [T]), we obtain 



2 2 



du 



dXr 



(A5) 



X 



a=l 13=1 

In the following, we drop the | attributes and keep in mind that there formally is a dependence 
on the position in the lattice. The expression for 5ai now becomes 

2 2 

niia 

(A6) 



1 ^ ^ 



a=l 8=1 



dXa 



Calculating explicit expressions for the quantities D^^o = di^dip leads to 



(A7) 



In a completely analogous way we obtain for 5c — recalling that c = ccz — the following 
expression (up to first order): 

duz 



be 



X + c + u{X + c) - (X + u{X)) 



c ^ c 



dz 



(A8) 



2. Corrections to the hamiltonian matrix 



The corrections to the hamiltonian involve the following quantities: 

^dV^ dV ^ 3a dV, , dV 

m.=\ ^^a, + —5c=- — {e.. + eyy) + c- 



1=1 

3 



da 

1=1 

Here, use has been made of the identity 

3 



3aM 1 



1=1 - V 1 

and the definition of the (first-order) components of the strain tensor: 

_ 1 / dui duj\ 



(A9a) 
(A9b) 



(AlO) 



(All) 



Note that the strain tensor is, in general, space-dependent: [e] = [e](X). Therefore, transforming 
the summation 



yddu^c^ = yci 



3adU dV 



(A12) 
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into a summation over g- vectors involves discrete Fourier transforms: 



^1 



q q' 



3adU , . ,_f dV ,_, . 

£xx{q -q)+ £yy{q -q)) + c-^e^^{q - q) 



2 da 



N ^ 

X. 



CAi(g ), 

(A13a) 
(A13b) 



Importantly, when the strain tensor [e] is space-independent, i.e. in the case of uniform deforma- 



tions, Eqs. (A13a) and (A13b) collapse into 



3a dU dU 



(A14) 



so that 



dU 
dc 



3a dU 

2 da di 
is g-independent. For the correction to 9)22 we then have 

3a dU 

09)22 = y^(exx + eTO). 



(A15) 



(A16) 



The off-diagonal non-zero hopping matrix elements are more complicated. We will need the 
quantities 



da ' a da 



1 dt ^ ^ duj3 



a=l 13=1 

dt I du^ 



ot± = -TT-oc = c 

dc dc dz 



(A17a) 
(Al7b) 



Let us consider the term 

3 



Xi 



1 1 dt 



2 2 



- /T7^5„Z]Z]I^I^^Ai®ea/3(9 -€)^a/3(9)cBi(g), (A18) 
1=1 ^ q q' a=l/3=l 



with 



Xi 



(Al9a) 
(A19b) 



du0 
dx 

1 dt 



Again, in the case of uniform deformations, is space-independent and Eq. (A18) simplifies to 

3 



a da 



q a=l f3=l 



(A20) 
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so that the matrix element 5i3i2 becomes 

a=l f3=l 



For F{q) one obtains 



^ ^ ^ ^ ^ 

„ .1 / icos(^ag^) + e-*i''''- ^isin(^ag„) , 
= 0^6*2'^^- 2 ^2 2 \ 2 ^yJ ^ _ ^^22) 

\ ^isin(^agy) |cos(^agj;) 

Near the K-point, one has 

F^^(g-i, + fc) = ^e^f I -|^ + 2(3^^^-M ^-|{k.-iky) ^ ^^^^^ 

The phase factor e*3 can be ehminated by redefining the creation and annihilation operators for 
electrons at B sites (position vectors Xi + di and X2 — di): 

ct - — ^cl -^e"'^, (A24a) 
cx.idl ^c^.id^'^- (A24b) 
Near the K-point, we obtain the following correction to the matrix element: 

6S)i2{k) = -i^^(exx[-l + '^{Sikx - ky)] + eyy[l + ^{ikx - 3ky)] + 2exy[i - ^{kx - iky) 

(A25) 

where the factor i comes from multiplying the phase factors e*3 and e*6 , and where the definition 



of the strain tensor in Eq. (All) has been used. For small deformations, the terms proportional to 



Ea/sk, with a,/3 G {x,y}, can be neglected: 



We now consider the term 

1 dt± 



E = E E 4. - «1^^^"'^A2 {q'), (A27a) 
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Here, the factor e*'' ''^ is equal to 1 because q' (only x- and y-components) and c (only z-component) 
are perpendicular. For space-independent we obtain 



^Jtxct c^^+.= c^^X]^Ai(9lcA2(g1 = c^e22X]cAi(9lcA2(g1, (A28) 




so that the correction to the matrix element i^is reads 




22 ■ 



(A29) 



These considerations then lead to the correction 5S) to the hamiltionan matrix S) given by Eqs. 
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